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Abstract 

In this paper, the space-time conservation ele- 
ment solution element (CE/SE) method [1] is 
tested in the classical axisymmetric jet insta- 
bility problem, rendering good agreement with 
the linear theory. The CE/SE method is then 
applied to numerical simulations of several 
nearly fully expanded axisymmetric jet flows 
and their noise fields and qualitative agreement 
with available experimental and theoretical re- 
sults is demonstrated. 

1 Introduction 

Jet noise is a challenging topic in computational aeroa- 
coustics. Generally, the computational scheme is re- 
quired on the one hand to resolve the acoustic waves 
without introducing too much dispersion and dissipation, 
while on the other hand, it is required to capture shocks, 
or other nonlinear phenomena, near or inside the jet cor- 
rectly. In addition, non-reflecting boundary conditions 
must be implemented in an efficient way either in the 
near field or in the far field. 

The recent ‘Space-Time Conservation Element and 
Solution Element Method' [2], or the CE/SE method in 
short, is a scheme that meets the above requirements. 
The CE/SE scheme possesses attractive properties for 
aeroacoustics computations in that: (i ) it possesses low 
dispersion and dissipation errors; (// ) its shock-capturing 
nature makes the computation of shock-cell structures 
simple and accurate; (///) the non-reflecting boundary 
conditions are simple and effective and can be applied 
in the near field of the jet without introducing excessive 
errors; and (iv) the vorticity, which plays an important 
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role in noise generating mechanisms, is determined with- 
out formal loss of accuracy. A detailed description of the 
CE/SE method can be found in the reports of Chang el 
al [3,4]. As demonstrated in our previous papers [5-8], 
the CE/SE scheme is well suited for computing waves in 
compressible shear flows [5], vorticity/shock interaction 
[ 5 - 7 ], as well as the near-field noise of an underexpanded 
axisymmetric supersonic jet w ith a shock cell structure 
[ 7 ] — all of the examples being comer stones of the jet- 
noise phenomena. 

In this paper, another aspect of jet aeroacoustics — the 
nearly fully expanded axisymmetric supersonic jet is in- 
vestigated numerically by using the CE/SE Euler solver. 
The paper is arranged as follows: The axisymmetric 
CE/SE scheme is briefly discussed in Section 2. Sec- 
tion 3 considers the classical jet shear layer instability 
problem and numerical results are compared to those ob- 
tained by linear normal-mode instability theory. In Sec- 
tion 4, the numerical results for several different types of 
noise fields of nearly fully expanded supersonic axisym- 
metric jets are presented and compared to experimental 
findings [10-12]. Conclusions are drawn in Section 5. 

2 Axisymmetric CE/SE Euler Solver 

In general, the CE/SE method systematically solves a 
set of integral equations derived directly from the phys- 
ical conservation laws. It is distinguished by the sim- 
plicity of its conceptual basis — flux conserv ation in both 
space and time. Because of its integral form of the phys- 
ical conservation laws, the scheme naturally captures 
shocks and other discontinuities in the flow 7 . Both de- 
pendent variables and their derivatives are solved for si- 
multaneously and, consequently, the flow vorticity can be 
obtained without reduction in accuracy. Non-reflective 
boundary conditions are also easily implemented be- 
cause of the flux-conservation formulation. 

2.1 Conservation Form of the Unsteady 
Axisymmetric Euler Equations 

Consider a dimensionless conservation form of the un- 
steady axisymmetric Euler equations of a perfect gas. Let 
p , u, i\ p , and 7 be the density, streamwise velocity com- 
ponent, radial velocity component, static pressure, and 
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constant specific heat ratio, respectively. The axisym- 
metric Euler equations then can be written in the follow- 
ing vector form: 

U t + FjC + Gy — Q, (1) 

where x, y ^ 0, and t are the stream wise and radial coor- 
dinates and time, respectively, and the conservative flow r 
variable vector U and the flux vectors in the stream wise 
and radial directions, F and G , are given by: 



where U m , F m , G m , m = 1 , 2 , 3, 4 are the same as in the 
standard 2-D CE/SE formulation [3, 4] and 



and U y are calculated at this node. Within a given so- 
lution element SE(j> n), where (j,fr,n) is the node 
index, the flow variables are not only considered con- 
tinuous but are also approximated by linear Taylor ex- 
pansions. This procedure is identical to the one in the 
standard 2-D CE/SE Euler solver [3, 4]. 

The discrete approximation of (2) is then 

<f h; ■ dS = V(CE) $UQm)l k , 

^ S ( CE ) (, j,k.n)eS(CE ) 

(3) 

for m — 1, 2, 3, 4, where = (F^, G* m , U* % ) denotes 
the linear Taylor-series approximation of the correspond- 
ing ‘ unstarred ' variables and Y^(jM,n)es(CE) = i * 
The right-hand side of (3), in general, is the volume of 
the CE times a weighted average of the ‘source’ term 
evaluated at the nodes on S(CE ). Each S(CE) is made 
up by surface segments belonging to two neighboring 
SE' s. All the unknowns are solved for based on these re- 
lations. No extrapolations (interpolations) across a sten- 
cil of cells are needed or allow ed. In order that the num- 


with 

Q i = —Us /y, Q2 = —U2U3/U1 y, 

Q 3 = -Ui/U lV , Q 4 = —Gi/y. 

The axisymmetric Euler equations can be thought of in 
this formulation as a variation of their two-dimensional 
counterparts with ‘source’ terms on their right-hand 
sides. Of course, these additional terms do depend im- 
plicitly on the solution to the equations and, hence, are 
not true source terms, but they can formally be treated as 
such for the purpose of deriving the numerical scheme. 

By considering (x, y , t) as coordinates of a three- 
dimensional Euclidean space F3 and using Gauss diver- 
gence theorem, it follows that Eq. (1) is equivalent to the 
following integral conservation law : 



where S( V) denotes the surface around a volume V in 
F3 and H m ~ (F m , G m , F m ). 

2.2 CE/SE Structure 

In the CE/SE scheme, the flux conservation relation in 
space-time is the only mechanism that transfers infor- 
mation between node points. The conservation element 
CE, or computational cell, is the finite volume to w hich 
the integral flux condition (2) is to be applied. Discon- 
tinuities are allowed to occur in the interior of a conser- 
vation element as long as the integrals on the right-hand 
side of (2) exist. A solution element SE associated with 
a grid node is here a set of three interface planes in F 3 
that passes through this node. The solution, i.e. U , U x , 
NASA/TM — 2000-2 1 0225 





Figure 1: Double staggered grid in F3. 


ber of equations derived from the flux conservation law 7 
above matches the number of unknow ns (here 1 2 scalar 
unknowns), the grid in (x, y) plane needs to be carefully 
designed [3, 4]. The mesh is staggered in both time and 
space. In a spatial plane in F 3 , the grid nodes, see Fig. 1 , 
are grouped as two staggered sets (open circles) and 
D2 (filled circles). At a given time step the unknowns 
are evaluated only on one of the grid sets, Q\ or and 
at the next it will be evaluated on the other set, i.e., the 
spatial sets alternate as time is stepped forward. As can 
be seen in this figure, each ‘interior’ node point then has 
three nearest neighbors at the previous time step. Thus, 
there are three CE ' s associated with each node point in 
this arrangement and, therefore, there are the same num- 
ber of relations as there are unknowns. 
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2.3 Treatment of the Source Term 

In solving (3), only the t 3™ k corresponding to the new 
(half) time level is taken to be nonzero (and equal to 
unity). This strategy has successfully been applied by 
Yu ef al [9] to problems with very stiff source terms. 

As the source term Q = Q(U) itself is a function 
of the unknown U at the new (half) time level, a local 
iterative procedure is needed to determine U. For our 
particular choice of 0j k , the discretized integral equation 
(3) reduces to the form 

U-Q(U)^=Uh, (4) 

where Uh is the local homogeneous solution (Q = 0 
locally). Note that Uh is in fact the standard 2-D Eu- 
ler solution at the previous (half) time step, obtained by 
explicit formulas. A Newton iterative procedure to deter- 
mine U is then 

jjd+D = V (i) _ _ u H ], 


the subject of many detailed analytical and experimental 
studies. 

The direct numerical computation of the shear-layer 
instability, using the CE/SE Euler scheme, is compared 
with linear results obtained using the normal mode ap- 
proach. The background mean flow r consists of a fast jet 
stream and a slow ambient stream around it, with the two 
parallel streams connected by a continuously changing 
annular shear layer. The streamwise flow variables in the 
fast jet stream are taken as the corresponding scales for 
the nondimensionalization. The length scale is taken as 
5/ 2, where S is the vorticity thickness and is defined as 

S — (T i* T'2*)/(dT */di/*) max 

with the subscript V here denoting dimensional quan- 
tities. In a parallel flow, the mean pressure is constant 
and the mean transverse velocity vanishes all over the 
domain, while the mean streamwise velocity and density 
profiles must be specified. 

In the linear stability theory, the normal-mode assump- 
tion is made, i.e.. 


where i is the number of iteration and 


4>(x, y, t ) = 6{y) exp[i(ax - art)], (5) 


*(£/) = U-Q(U) — . 


Normally, U at the previous (full) time step is a good 
initial guess U 1 ^ 0 - and the procedure takes about 2-3 it- 
erations to convergence. The Jacobian matrix is given 
by 


m 
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The inverse of the Jacobian, i.e., (f§) _1 can easily be 
derived analytically for this particular case, thus, leading 
to savings in CPU time. 


3 Free Jet Shear Layer Instability 

In the first example, we study the inviscid linear and 
nonlinear instability properties of an axisymmetric free 
shear layer (jet). This important class of flows forms 
the basis of jet noise generation theory and have been 


where 4>{x, y, t) denotes the perturbation streamwise and 
radial velocities u\ v\ pressure p\ or density p f \ and a, 
w’ are the complex wave number and real angular fre- 
quency of the disturbance, respectively. We note that u, 
t\ p and p are also complex functions in general. Adding 
the above perturbations to the corresponding mean flow 
variables and substituting them into the 3-D axisymmet- 
ric Euler equations, a set of equations for the unknown 
perturbation flow variables is obtained. Assuming that 
the mean flow variables satisfy the Euler equations and 
by eliminating those terms of second order or higher, the 
compressible Rayleigh equation in terms of the perturba- 
tion pressure eigenfunction component p can be achieved 
after some manipulation. 

To be specific, the streamwise mean flow ? velocity U is 
here taken as 


U(y) = 


1 + R tanh y 
1 -f- R 


( 6 ) 


where R is the velocity ratio of the shear layer. We 
further assume for simplicity that the Prandtl number is 
unity and the mean flow temperature T distribution is 
then obtained from the Buseman-Crocco relation: 

T = T 2 +^l(U-U 2 )+l('r-l)Mf(U-U 2 )(l-U), 

I — 02 

( 7 ) 

where the subscripts 1 and 2 denote the fast stream and 
slow stream variables respectively. T in turn yields the 
density distribution across the shear layer. As in [6], we 
choose the parameters 


R = 0.15, J\ij = 1.5, 72 = 1.85 
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Then, the nondimensional flow states are worked out as: 
V\ = 1 , V\ = 0 , Px = 1 / 3 . 15 , 

Tx = 1, pi = 1, Mi = 1.5, 

[/a = .7391304, V 2 = 0, p 2 = 1/3.15, 

T 2 = 1.85, p 2 = 0.5405405. 

Under the above conditions, with an appropriate a cho- 
sen, the complex eigenfunction component p = p r 4- ipj 
and the corresponding eigenvalue a = a r -f iaj, where 
the subscripts r and i denote real and imaginary parts 
respectively, can be found by solving the compressible 
Rayleigh equation. 

3,1 Linear, Nonlinear Instability of an 
Axisymmetric Jet and Vortex Roll-Up 

In studying the instability of a free jet shear layer and the 
consequent vortex roll-up and noise generation, one is 
mainly interested in the most unstable mode which dom- 
inates the nonlinear vortex roll-up. In order to investi- 
gate such instabilities, we shall enforce a small harmonic 
perturbation at a single frequency at the inlet boundary 
and compute the resulting development of the spatial in- 
stability wave. In this test example, the computational 
domain spans between 0 ^ x ^ 300 and 0 ^ y ^ 20, 
with a grid of 601 x 101, respectively (ie. Ax = 0.5 and 
Ay = 0,2 ). A time step size At = 0.15, e = 0.5 and the 
weighted average index a = 0 were found appropriate. 
The computation of the unsteady jet flow is carried out 
until t = 390 ensuring that the spatial instability is fully 
developed. 

From a chart of a versus -a*, for the case under con- 
sideration here, it can be observed that the exponential 
growth rare - a* reaches its maximum of about 0.0244 
when the angular frequency w = 0.405 and the latter is 
the value used in the following test example. The com- 
plex spatial eigenvalue a = 0.44726 - 0.02444i and 
the corresponding eigenmode is then determined by solv- 
ing the compressible Rayleigh equation for this value of 
u). The eigenfunction is normalized such that the max- 
imum r.m.s value of the streamwise velocity component 
is unity. The forcing at the inlet boundary (x = 0) then 
has the form 


<f> = 6[4> r cos art + 4>i smart], 


( 8 ) 


d<p 

dx 


6[—(<t> r ai + 4>%<Xr ) cos art — (&c*i ~ <£ r a r ) sin art], 


dy 


— <5 [4> r cosurt -h $ smart), 


(9) 

(10) 


demonstrate the capability of the CE/SE scheme to com- 
pute both linear and nonlinear instability waves, the forc- 
ing amplitude 5 is taken to be 1/1000. 

Fig. 2 displays the flow variables contours. Vor- 
tex roll-up is clearly observed in this figure. This fig- 
ure clearly demonstrates the effectiveness of the non- 
reflecting boundary conditions at the top, bottom and 
outlet. We emphasize that the domain shown in the 
figures is exactly the computational domain, no buffer 
zones, cut-offs or other numerical fixes are applied. 

To assess the grid dependence of the numerical results, 
the computation was repeated with half grid sizes and 
time step size. Fig. 3 shows that they are quite similar 
except that, as expected, the finer grid allows higher har- 
monics to he better resolved. 

Fig. 4 shows the power spectral density PSD of the 
computed p f in log 10 scale at the streamwise station x = 
150. The spectral peak at the forcing frequency / = 
~ = 0.06446, the second, third, and fourth harmonic 
thereof are clearly displayed. 

c«nvpfd 



Figure 4: Power spectral density of the jet shear layer at 
x = 150. 


4 Mach radiation of jets 

In Oerters experiments [10-12] and Tam and Hu’s 
analysis [13], three distinct types of Mach waves for 
high-speed jets were identified. Our purpose in this sec- 
tion is to investigate these three types of axisymmet- 
ric fully expanded circular jets. The numerical results 
are presented and qualitatively compared to experimen- 
tal findings [10-12]. In all the examples, a perturbation 
(or forcing) is provided by a right-hand-side source term 
in the energy equation, the 4th component of Eq. ( I ), lo- 
cated at the origin (0, 0) inside the jet core: 


where <f> denotes any of the eigenfunction components 
u\ v\ p\ or p\ and <5 is the forcing amplitude. To fully 


^ exp[-J5(x 2 + y 2 )] cos (art) 
7-1 
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where 6 is a small number (0.00005 < 6 < 0.001), 
uj = 2 nSt, and the constant B = 8. Even though 
the CE/SE implementation of non-reflective boundary 
conditions is very effective (generally less than 1% re- 
flection), it is of advantage to also add a (streamwise) 
buffer zone to the outflow' side of the computational do- 
main to further reduce numerical reflections from radial 
regions where the mean flow 7 changes rapidly, A buffer 
zone of 20 cells with exponentially increasing size was 
used in the computations presented below and effectively 
kept the computed acoustic field free of such reflections. 
Due to the intrinsic properties of the CE/SE scheme, it 
is likely that the number of buffer region cells could be 
substantially reduced without a significant degradation of 
the computed results — this is yet to be studied, however. 

4.1 Mach radiation from a supersonic 
axisymmetric jet 

In this test example, a supersonic jet with Mach num- 
ber Mj = 2.0 is considered. As the jet flow is fully 
expanded, the jet nozzle exit pressure is the same as the 
ambient pressure. However, tiny fluctuations of pressure 
are still allowed and accounts for the acoustic w 7 ave gen- 
eration and propagation. The overall motion is approx- 
imately in an axisymmetric mode [14] at the Strouhal 
number St = 0.2. We note that since the Eiiler solu- 
tion does not account for the shear layer spreading due 
to viscous effects, our numerical results are not expected 
to perfectly match the experimental [15] sound pressure 
level, SPL. 

The computational domain spans between 0 ^ x ^ 
33 D and 0 ^ y < 19D, with a non-uniform grid of 
350 x 281 nodes, where D is the diameter of the jet at 
the nozzle exit. More grid points are packed around the 
shear layer. The last 50 cells in x direction grow ex- 
ponentially in length and are used as buffer zone. The 
perturbation source strength 5 = 0.001. A time step size 
At — 0.005, e = 0.5 and the weighted average index 
a = 0 are found appropriate. The computation of the 
unsteady jet flow is carried out until t = 100 ensuring 
that the spatial instability is fully developed. 

Fig. 5 illustrates the isobars and v-velocity contours 
in the near and intermediate field. The Ma.ch_ radiation 
flow pattern agrees qualitatively with experimental [15] 
and other computational [14] results and is categorized 
by Oertel [10] as Type I Mach radiation. 

4.2 Mach radiation from a supersonic 
axisymmetric jet at low Strouhal number 

In this test, the same supersonic Mj = 2.0 jet is consid- 
ered but a low Strouhal number St = 0.07 is chosen. The 
computational domain spans between 0 ^ x ^ 8D and 
0 < y ^ 2.5D, with a non-uniform grid of 630 x 151 
nodes. The last 30 cells in x direction grow 7 exponen- 
tially in length and are used as buffer zone. The per- 


turbation source strength 5 = 0.0001. A time step size 
At = 0.002, e = 0.5 and the weighted average index 
a = 0 are found appropriate. The computation of the 
unsteady jet flow is carried out until t = 27. 

Fig. 6 demonstrates the isobars in the near field. At a 
low Strouhal number, the radiated Mach weaves are quite 
weak and rather steeply inclined. They look almost ver- 
tical. This is categorized as Type II Mach radiation by 
Oertel [11]. Tam and Hu [13] also give an explanation 
of the phenomena in their analysis. Inside the jet core, 
the cross-hatch pattern of a trapped Mach wave system 
is observed. As will be seen in the next subsection, such 
a phenomenon is more pronounced when Mj is low . 

4.3 Mach wave system inside the jet 

This is the Type III experimental result by Oertel [12]. 
It was found that, as analyzed by Tam and Hu [13], the 
Mach radiation attenuates exponentially outside the jet 
shear layer, while the Mach wave system is trapped in- 
side the jet and forms a cross-hatch pattern. Here the jet 
shear layer plays a role as a ‘w 7 ave guide'. We consider 
two cases in this example. 

4.3.1 Case 1: Mj = 1.2 In the first case, a su- 
personic jet with low Mach number Mj = 1.2 is con- 
sidered. The computational domain is 0 < x < 3D and 
0 < y ^ ID, with a non-uniform grid of 330 x 101 
nodes, where again D is the diameter of the jet at the 
nozzle exit. The last 30 cells in x direction are used as 
buffer zone. Here, the near field is chosen for investi- 
gation since our attention is focussed on the Mach wave 
system inside the jet. The perturbation source strength 
5 = 0.00005, with a high Strouhal number St — 2.0. A 
time step size At = 0.0015, e = 0.5 and the weighted 
average index a = 0 are chosen. The computation of the 
unsteady jet flow 7 is carried out for 1 2000 steps to ensure 
that the spatial instability is fully developed. 

Fig 7 clearly demonstrates the cross-hatch pattern of 
the Mach wave system inside the ‘wave guide’ — the jet 
core. The Mach waves are repeatedly reflected because 
they do not penetrate the annular shear layer. The pattern 
qualitatively agrees with Oertel's experiment [12] (Type 
III) and the analysis of Tam and Hu [13]. At a low su- 
personic Mach number Mj = 1.2, at the outer surface 
of the shear layer, the local phase velocity can hardly ex- 
ceed the local speed of sound and hence practically no 
exterior Mach radiation occurs. 

4.3.2 Case 2: Mj = 2.0 In the second case, the 
Mach number Mj is increased to 2.0. The computational 
domain is 0 ^ x 7.5D and 0 ^ y ^ ID, with a non- 
uniform grid of 280 x 101 nodes. As before, the last 
30 cells in x direction are used as buffer zone. The per- 
turbation source strength 6 = 0.001, with the Strouhal 
number St = 2.0. A time step size At = 0.002, e = 0.5 
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and the weighted average index a = 0 are chosen. The 
computation of the unsteady jet flow is again carried out 
for 12000 steps. 

Fig 8 displays the cross-hatch pattern of the Mach 
wave system inside the ‘wave-guide’ jet core. The pat- 
tern qualitatively agrees with the analysis of Tam and Hu 

[13]. At a supersonic Mach number Mj = 2, the local 
phase velocity at the outer shear layer may exceed the 
local speed of sound and hence Mach radiation is gener- 
ated. 

5 Concluding Remarks 

In this paper, the noise field of nearly fully expanded 
axisymmetric jets are investigated by using the recent 
CE/SE Euler scheme. For the classical jet-shear-layer 
instability problem, the numerical results show a linear 
initial development of the spatially growing instability 
wave (in good agreement with stability theory 7 ) followed 
by nonlinear effects leading to vortex roll up. Many as- 
pects of the computed results for Mach radiation are also 
in good agreement with Oertel's experimental findings 
[10-12] and Tam’s analysis [13]. 

The advantages of the scheme as claimed in Section 
1 are confirmed in the present simulations. The imple- 
mentation is ‘effortless’ in that no special treatment and 
parameter selections are needed. 
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Figure 3: Comparison of density contours for coarse and fine grids at t = 390 
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Figure 5: Type I Mach radiation of a supersonic jet at Mj = 2,0 and Strouhal number St = 0.2, domain size 33Dxl9D, 
non-uniform grid. 
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Figure 6: Type II Mach radiation of a supersonic Mj = 2,0 jet at low Strouhal number St=0.07, domain size 8Dx5D, 
non-uniform grid 630x150. 



MJ=1.2, 3Dx1D, 300x100 grid, St. =2. dt=.0015,ee=.00005; 
showing cross-hatching Mach wave system. 

Figure 7: Mach wave system within a supersonic jet of Mj = L2 at high Strouhal number St = 2.0, isobars show the 
cross-hatch pattern. 
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Mach radiation periodic Mach wave system 



Figure 8: Type ITT Mach wave system within a supersonic jet of Mj = 2 at high Strouhal number St 
show also weak Mach radiation. 
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